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ABSTRACT: This paper presents an original ABC algorithm, ABC 
Shadow, that can be applied to sample posterior densities that are con¬ 
tinuously differentiable. The proposed method uses the ideas given by 
the auxiliary variable MH of [18]. The obtained algorithm solves the 
main condition to be fulfilled by any ABC algorithm, in order to be use¬ 
ful in practice. This condition requires enough samples in the parame¬ 
ter space region, induced by the observed statistics [8]. The algorithm is 
tuned on the posterior of a Gaussian model which is entirely known, and 
then it is applied for the statistical analysis of several spatial patterns. 
These patterns are issued or assumed to be outcomes of point processes. 
The considered models are: Strauss, Candy and area-interaction. 
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1 Introduction 

Let us assume that some spatial data set is observed. Typical examples of 
such sets are digital images, epidemiological or environmental data or cata¬ 
logues of celestial bodies in astronomy. One typical question related to these 
examples is the detection and the characterization of the “hidden” pattern 
in the data. Such patterns may be the collection of cells in some biological 
image, the set of clusters exhibited by a surveyed disease or the filamentary 
network outlined by the galaxy positions within the observed Universe. 

Within a probabilistic context, this problem is tackled by assuming that 
the pattern is the outcome y of a stochastic process Y. A possible solu¬ 
tion within this context is probabilistic modelling and maximisation. This 
option requires as a second hypothesis, the knowledge of the model param¬ 
eters. The model choices may include among others, random fields [25] or 
point processes nailg]. 

Let us consider the observation window to be a finite domain VU C with 
volume 0 < v{W) < oo, where v is the Lebesgue measure on The hidden 
pattern is supposed to be the realisation y of a marked point process Y on 
W X M. This means that the hidden pattern is a random set made of a 
finite number k of random objects y = {yi,y 2 , • • • j 2 /k}- These objects have 
a location parameter in W and a characteristic vector or a mark in M, that 
is yi = G VL x M for z = 1,..., k. The marks probability space is 

{M,}A,vm)- The pattern probability space is /r) with 0 the configu¬ 

ration space, T the associated it— algebra and y the reference measure. 

The reference measure y, is given by the unit intensity marked Poisson point 
process. This process is defined as follows: the number of objects is chosen 
with respect to a Poisson law of parameter z^(LP), then the objects loca¬ 
tions Wi are independently and uniformly chosen in W, and finally to each 
location a mark rrii is attached. The marks are independent of the objects 
locations and independent of the other marks, and they all follow the same 
distribution vm- 

Due to the independence properties of the process, there is no interaction 
among the objects in a configuration. Therefore, the outcome of such a 
process is considered to be a completely random configuration of objects, a 
completely random pattern. More structured patterns may be obtained by 
introducing interactions among objects. Such interactions may be specified 
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by means of a probability density p with respect to the reference measure 
/i. The Gibbsian modelling framework allows to write such a probability 
density as 


p{y\^) 


exp[-?7(y|6*)] 

0 ( 9 ) 


( 1 ) 


with [/ : M'*' the energy function, 9 the model parameters and c{9) the 

normalising constant. 


Being in the possession of a model ([T]), the hidden pattern estimator is given 

by 

y = arg max{p(y|0)} = arg min{C/(y|0)}. (2) 

yeo yen 

The optimisation in ([2]) is done using a simulated annealing algorithm based 
on a Metropolis-Hastings (MH) or spatial birth-and-death dynamics [111121] . 
It is important to notice that these algorithms do not require the compu¬ 
tation of the normalising constant c{9) which is not always available in 
analytical closed form. 


The dual formulation of the pattern detection question is the parameter 
estimation problem. Let us now consider that an object pattern y is observed 
in W. The observed pattern is supposed to be the realisation of a marked 
point process given by the probability density p{y\9). Let p{9\y) be the 
conditional distribution of the model parameters or the posterior law 

pwy) = ( 3 , 

p{9) the prior density for the model parameters and Z{y) the normalising 
constant. The posterior law is defined on the parameter space 0. For sim¬ 
plicity, the parameter space is considered to be a compact region in with 
r the size of the parameter vector. The parameter space is endowed with its 
Borel algebra 7e. 

Straightforward sampling of ([3]), for instance using a classical MH dynam¬ 
ics, cannot be always done, since it may require the evaluation of the ratio 
c{9)/c{'ip) for (0,'!/’) £ 0^- The theoretical solution to this problem is given 
by [18]. The authors propose an auxiliary variable MH sampler that, by a 
clever and elegant choice of the proposal and the auxiliary variable densi¬ 
ties, avoids the computation of the normalising constants. Nevertheless, as 
the authors themselves explain in detail, it is difficult to show how these 
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choices prevent the simulated chain from poor mixing. Approximate solu¬ 
tions for posterior sampling are provided by Approximate Bayesian Compu¬ 
tation (ABC) methods O [9l [151 [1] . These methods are very attractive from 
a practical point of view. The theoretical developments of ISE] established 
nice connections of the ABC based inference with kernel and /c—nearest 
neighbour estimation. The bias and the variance for the posterior distri¬ 
bution estimate is given by [8]. Despite its simplicity, the main criticism 
against ABC methods is the large number of generated samples that cannot 
be used. In fact, in order to get reliable numerical results, an ABC method 
should produce enough samples in the parametric region induced by the 
posterior sampling. 

The aim of this paper is to solve this drawback for posterior densities that 
are continuously differentiable. This a strong but realistic assumption. With 
this goal in mind, we embed the auxiliary sampling ideas of [18] in the ABC 
framework given by The paper continues with the presentation of the 

auxiliary variable method of [T8|. Next the ABC methods from I31I1 are 
described. The principles of our posterior sampling method are shown in 
Section 3. The application results are analysed in Section 4. The method is 
first tuned on a Gaussian model for which the posterior density is entirely 
known. Then the proposed algorithm is applied to the statistical analysis 
of spatial patterns issued or assumed to be point processes outcomes. The 
considered models are: the Strauss, the Candy and the area-interaction pro¬ 
cesses. At our best knowledge, it is the hrst time an ABC based statistical 
analysis is done for such models. The paper ends with conclusions and per¬ 
spectives. 


2 Sampling posterior distributions 

2.1 Theoretical solution: auxiliary variable MH algorithm 

The goal is to sample the posterior 

p{e\y) (xp{y\e)p{e). 

The authors idea in |18j is to use an auxiliary variable X with probability 
density a(x|0,y). Consequently, the joint distribution becomes 

^(6',x|y) = o(x|6l,y)p(6'|y) = q(x| 6I, y) 


(4) 
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Sampling from ([5]) can be done using a MH algorithm that uses a proposal 
made of two components, such as: 

qii9,x) (6l',x')) = qi{e'\9,x)q2{x'\e',9,x). 


The first component qi may be a very simple probability law and it may not 
depend on x. Let us consider, for instance, the symmetric distribution 

q,i9'\9)=qi{\9'-9\). (5) 


The second component q 2 is chosen independently on (0,x) as follows 

exp[—C/(x'|0')] 


g 2 (x \9 ,6l,x) = g 2 (x \9 ) = 


c{9') 


( 6 ) 


The MH ratio H£){{9, x) —)■ {9', x')) of the induced dynamics is, by definition, 


i/D((0,x)^(0',x')) 


7 ^(g^xV) g((g^xO^(g,x)) 

7r(6',x|y) q{{9,x) ^ {9',x'))' 


(7) 


Plugging the joint density ([4]) and the proposals ([5]) and ([6]) into the 
MH ratio becomes 

ai^'\9\y)eM-Uiy\0'M0’) exp[-f/(x|g)] 

a{x\9, y) exp[-f/(y |6»)]p(6') exp[-U (x'|6»')] ’ 

which does not depend on the normalising constant c. 


For exponential family models, the MH ratio has a simpler form. In this 
case, the energy functions are of the form U{y\9) = {t{y),9), where (,) is 
the scalar product of the sufficient statistics vector of the model t{y) and 
the parameter vector 9, so the MH ratio is 

a{yL'\9', y)p{9') exp[-(t(y), 6*^] exp[-(t(x), 6>)] 
a{x\9, y)p{9) exp[-(t(y), 9)] exp[-(t(x'), 6»')] 

Under these considerations, the auxiliary variable MH sampler works as 
follows: 


Algorithm 1. Assume the observed pattern is y and the current state is 
{9kj Xfc). 


1 . 

2 . 


Generate a new candidate 9' using the proposal qi{9'\9). 


Generate a new candidate x' using the proposal q2{'x'\9') 


exp[—C/(x'| 0 ')] 

W) 
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3. Compute the MH ratio H£){{6,x) —>■ (0',x')) using ([7]) or one of its 
simplified forms previously detailed. 

4- The new state (0^+1) = (0^x') is accepted with probability 

a = min{l,i7^((6',x) (^^x'))} , 

otherwise, the chain remains in the same initial state, that is 

(^A:+l)^fc+l) — {dkjXjf,'). 

5. If more samples are needed, re-iterate the whole procedure. 

The outputs of the Algorithm [1] converge in law towards the unique equi¬ 
librium distribution 7r(0,x|y) [18]. This result holds if sampling from the 
proposal q 2 {x'\ 6 ') is done exactly. Despite the fact that the convergence time 
of the exact simulation methods for marked point processes may highly de¬ 
pend on model parameters [Tl|, we do not consider this characteristic as a 
drawback of the auxiliary variable method, since this method is mathemat¬ 
ically correct. We hope that future advances in computer science and algo- 
rithmics will be able to propose a practical implementation of this solution. 
The difficulty in a straightforward use of the auxiliary variable method is 
the freedom of choice of the auxiliary variable density. As the authors them¬ 
selves showed in [T8|, there is no clear choice for it, such that good mixing 
properties of the simulated chain are guaranteed. A possible strategy, that 
may be fair enough, is to set a{x\0,y) = where 6(y) is the 

pseudo-likelihood estimate of the model parameters based on the observa¬ 
tion of y. 


2.2 Approximate solution: ABC methods 

ABC (Approximate Bayesian Computation) is the generic name for numer¬ 
ical simulation methods allowing statistical inference based on an approxi¬ 
mate sampling from the posterior distribution p(0|y) [ll[5l[9l[T5]. Its general 
idea is described by the following algorithm. 


Algorithm 2. Assume the observed pattern is y, fix a tolerance threshold 
e and an integer value n. 


1. For i = I to n do 
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• Generate 6 i according top{9). 

• Generate Xj according to the probability density 

p(x|ft) = 

^\9i) 

2. Return all the 9i’s such that the distance between the statistics of the 
observation and those of the simulated pattern is small, that is 

d{t{y),t{^i)) < e 

In the general case, the outputs of the Algorithm [2] are distributed accord¬ 
ing to the law 7r(0|d(t(y), t(x)) < e). The author in [8] established a link 
between kernel estimation and inference based on this algorithm. In this 
paper, estimates of the bias and the variance of the posterior density esti¬ 
mate are given. The choices influencing the algorithm outputs are related 
to the statistics vector, the distance definition and the tolerance threshold. 
For exponential family models the sufficient statistics vector is the natural 
choice. 

The authors in [7j considered a slightly different version of the previous 
algorithm. 

Algorithm 3. Assume the observed pattern is y and two integer values kn 
and n such that 1 < kn < n. 

1. For i = 1 to n do 

• Generate 9i according top{9). 

• Generate Xj according to the probability density 

c{9i) 

2. Return all the 9i’s such that t{xi) is among the kn-nearest neighbours 
oft{y). 

The Algorithms [3] and [2] are similar. The same requirements hold for both 
algorithms. But for the nearest neighbour algorithm no tolerance threshold 
is needed. Instead, the kn parameter should be fixed. Intuitively, the smaller 
e is or the higher kn is (such that kn/n 0 with n oo), the closer the 
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algorithms outputs are to the posterior density. In both cases, the posterior 
density is estimated using a kernel of a given bandwidth, over the selected 
samples. Clearly, the bandwidth of the kernel influences the quality of the 
estimation. This is not a drawback in our opinion, whenever enough points 
are sampled close to the observed statistics. So, as indicated in [8], the 
key point in building an ABC procedure is to ensure that the sampling 
mechanisms of the algorithm do not produce too many samples that are 
sparse around the observed statistics. 

3 ABC Shadow algorithm 

The ABC method we propose follows directly the recommendation of [8], 
to sample close to the posterior. Our solution is inspired by the auxiliary 
variable method [18]. This goal is achieved by building two Markov chains, 
the ideal and the shadow chains. The ideal chain is a theoretical chain that 
cannot be used in practice, but its equilibrium distribution is the posterior 
law we want to sample from. The shadow is a chain that can be practically 
simulated, and it follows as close as desired the ideal chain during a fixed 
number of steps. 

3.1 Ideal chain 

In theory, Markov chain Monte Carlo algorithms may be used for sampling 
p{0\y). For instance, let us consider the general MH algorithm. Assuming 
the system is in the state 9, this algorithm first chooses a new value V’ 
according to a proposal density q{6 —)■ tp). The value ip is then accepted 
with probability ai{9 —>■ ip) given by 




The transition kernel of the Markov chain simulated by this algorithm is 



with A G 7e. 
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The conditions that the proposal density q{9 —>• ijj) has to meet, so that the 
simulated Markov chain has a unique equilibrium distribution 

tt{A) = [ p{e\y)dij{e), 

Ja 

are rather mild [23]. Furthermore, the simulated chain is uniformly ergodic, 
that is, there is a positive constant M and a positive constant p < 1 such 
that 

sup II -vr(-) ||< Mp*", neN. 

0e0 

For a hxed A > 0, a parameter value € 0 and a realisation x of the model 
p{-\v) given by ([T|), let us consider the proposal density 

q{e -aiIj) = q/:^{e -A iplx) = ( 9 ) 

with /(xlV’) = exp[—f7(x|'!/))]. Here lb( 6 i,A/ 2 ){'} is the indicator function over 
6(0, A/2), which is the ball of centre 0 and radius A/2. Finally, 7(0, A, x) 
is the quantity given by the integral 

I{e,A,x)=f f{yL\(j))/c{(j))d(p. 

Jb{e,A/2) 

This choice for q{9 -A ip) guarantees the convergence of the chain towards vr 
and avoids the evaluation of the normalising constant ratio c{9)/c{ip) in ([8]). 
We call the chain induced by these proposals the ideal chain. Nevertheless, 
the proposal Q requires the computation of integrals such as 7(0, A, x), 
and this is as difficult as the computation of the normalising constant ratio. 
Still, this construction allows a natural approximation of the ideal chain: 
the shadow chain. 


3.2 Shadow chain 

In order to ease the reading of this section, the proofs of the theorems and 
propositions are given in the Appendix at the end of the paper. 

Let Va denote the volume of the ball 6(0, A/2), and let 

= ^l6(0,A/2){V'} 


Ua{ 0 -0 iP) 


( 10 ) 
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be the uniform probability density over the ball b{9, A/2). The probability 
density p(x|</>) defined by ([T]) is assumed to be a continuously differentiable 
function in 0 , p(x|-) G C^( 0 ). 


Theorem 1. Let x be a point in such that the function p{'x\(l)) is strictly 
positive and continuous with respect to cf, then we have that: 

(i) The probability distributions given by the proposal densities qA{9 ■) 
and Ua{9 —> •) are uniformly as close as desired in 9 as A approaches 
0. That is, for any fixed 9 € Q and A £ Te, we have 

lim [ \qAi9 i/) — Ua{9 —)• = 0 . 

^-^ 0 + Ja 


(a) For any fixed 9 £ &, the quotient functions and 


^h(e,A/ 2 )(') 


are uniformly as close as desired in 9 as A approaches 0. That is, for 
any fixed 9 £ &, we have 


lim sup 

^->■ 0 + i/)g 0 


qA{9 -A V^lx) 
qAitp ^'Ix) 


^lK»,A/ 2 )W) 


uniformly in 9 £ Q. 

Moreover, i/p(x|-) G C^(0), then rates of convergence in (i) and (ii) can be 
provided. 


This result leads to the construction of a new Markov chain that approxi¬ 
mates the behaviour of the ideal chain for small values of A. The first part 
of the theorem allows the use of the uniform law Ua{9 -a if) instead of 
qA{9 -A if) for proposing new values. The second part of the result enables 
the approximation of the ideal acceptance ratio. The resulting new chain 
is called the shadow chain. Assuming the chain in a state 9, a new value 
if is chosen uniformly in the ball 6(0, A/ 2 ). The state if is accepted with 
probability 


as{9 -A if) = min 


. Pj'ifly) /(x|0)c(V>)lfe(.0,A/2){^} ) 

’ P(6'|y) ^ /(x|V')c(0)lfe(0,A/2){V'} J ’ 


( 11 ) 


otherwise the chain remains in the initial state. 
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By construction, the shadow chain is irreducible and aperiodic ([IZ]). Yet, 
we do not have any knowledge about the existence and the uniqueness of its 
equilibrium distribution. The following corollary is a direct consequence of 
Theorem [TJ 


Corollary 1. The acceptance probabilities of the ideal and shadow chains 
given by ([8]) and dm) respectively, are uniformly as closed as desired when¬ 
ever A approaches 0. 

Next, we show that, for a given x, it is possible to choose A such that during 
n steps, the ideal and shadow chains evolve as close as desired. 

Proposition 1. Let Pi and Pg be the transition kernels for the ideal and 
the shadow Markov chains using a general A > 0 and a configuration x G 
as in Theorem [IJ respectively. Then for every e > 0 and every n G N, there 
exists Aq = Ao(e, n) > 0 such that for every A < Aq we have — 

Ps'^\6,A)\ < e uniformly in 6 £ Q and A G 7e. If p{'x.\9) G C^(0), then a 
description o/Ao(e, n) can be provided. 

The previous results allow the construction of the following ABC algorithm. 
Its outputs are approximate samples from the posterior p{6\y). 


Algorithm 4. ABC Shadow : fix A and n. Assume the observed pattern 
is y and the current state is 9 q . 


1. Generate x according to p{x.\9o). 

2. For k = 1 to n do 


• Generate a new candidate -ip following V’) defined 

by (Ho]). 

• The new state 9k = ip is accepted with probability as{9k-i —>■ ip) 
given by (HU), otherwise 9k = 9k-i. 

3. Return 9n 


4 . If another sample is needed, go to step 1 with 9 q = 9n. 
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If n ^ cx) the Algorithm 0] may not follow closely an ideal chain started in 
Oq. In this case, the shadow chain drifts away from the stationary regime 
of the ideal chain. Nevertheless, for a fixed n and A —>■ 0, the ideal chain 
approaches the equilibrium regime, equally fast from any initial point. If e, 
the distance after n steps between the ideal and the shadow chain is small 
enough, the shadow chain will also get closer to the equilibrium regime. The 
triangle inequality gives a bound for the distance after n steps, between the 
shadow transition kernel and the equilibrium regime 

||pW(0,.)-vr(.)|| <M(x,AK + e. 

One very important element here is the auxiliary variable x. Refreshing it 
allows to re-start the algorithm for another n steps more, and by this, to ob¬ 
tain new samples of the approximate distribution. Comparing with classical 
ABC, no useless samples are produced, in the sense that the outputs of the 
Shadow algorithm tend to have a distribution “closer” to the true posterior 
while X is renewed. In the next section, a simulation study is done, in order 
to analyse the Shadow algorithm performances. 


4 Applications 

4.1 Method analysis: posterior approximation for a Gausian 
model 

The Algorithm 0] can be applied to exponential family models. Here, the 
method is tested on the posterior of a Gaussian model. In this case, the 
posterior distribution has an analytical expression allowing its simulation 
using a MH dynamics. The comparison of the ABC Shadow algorithm with 
the classical MH dynamics allows the tuning of the algorithm and its numer¬ 
ical analysis in terms of mixing, parameter sensibility and initial conditions 
dependence. 


The posterior distribution of a Gaussian model with mean ffi and variance 
02 is 


p(0i,02\y) oc 


exp (^G(y) - 

0 ( 01 , 02 ) 



Pi01,02), 


( 12 ) 


with t(y) = (ti(y), t 2 (y)) the observed sufficient statistics vector. If a sample 
of size m is observed, that is the sequence y = {Yi(a;), ¥ 2 ( 0 ;),... , Ym(w)} is 
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given by m independent random variables following a Gaussian distribution 
with parameters 9i and 62 , then the sufficient statistics vector is 

( m m 

i=l i=l 

For our experiment, we simulated m = 1000 independent identical Normal 
random variables with parameters 9 = = (2,9). The observed suf¬ 

ficient statistics vector is t{y) = (1765.45,12145.83) and p{9) the uniform 
distribution over the interval [—100,100] x [0,200]. The MH algorithm pa¬ 
rameters implemented to sample from (1121) used uniform proposals of width 
(0.5,0.5) around the current values. Finally, this dynamics was run for 
125 X 10^ iterations and samples were kept every 12500 steps. This gave an 
amount of 1000 samples. For the ABC Shadow algorithm, the A parameter 
was set to (0.005,0.025) and n = 500. The algorithm was run 25 x 10^ times. 
Samples were kept every 25 repetitions of the ABC procedure. This gave an 
amount of 1000 samples as for the MH dynamics. 

Figure [T] presents the results of the MH and the ABC Shadow algorithms 
used to sample from the Normal posterior (11211 . Table [T] shows the values 
of some summary statistics obtained using both methods. The posterior 
approximation looks satisfactory. It is important to notice, that the mean 
and the median posterior estimates are close to the maximum likelihood 
estimates of the model parameters. 



Summary statistics for Normal posterior sampling 

Algorithm 

Qb 

Q 25 

QbO 

0 

Q 75 

Q 95 

MH 9i 

1.60 

1.69 

1.75 

1.76 

1.82 

1.92 

ABC 9i 

1.60 

1.70 

1.76 

1.76 

1.82 

1.91 

MH 02 

8.45 

8.80 

9.07 

9.08 

9.33 

9.76 

ABC 02 

8.35 

8.78 

9.03 

9.06 

9.33 

9.83 


Table 1: Empirical quantiles and mean for the posterior of the Normal 
model. 

Figure [2] shows the sampling paths obtained using the MH and the ABC 
Shadow algorithms, respectively. In both cases, the sample paths exhibit a 
similar behaviour. The ABC chain indicates rather good mixing properties, 
in the sense that the chain does not get locked in a certain parameter space 
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Normal posteriors - second parameter 


Normal posteriors - second parameter 



Figure 1: Normal posterior sampling using a MH algorithm and the ABC 
Shadow procedure; boxplots and qqplots of the algorithms outputs. 
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region and that it covers the parameter region given by the posterior distri¬ 
bution. 

Summary statistics obtained for the ABC algorithm run for different pa¬ 
rameters A and initial conditions 0q are presented in Table [21 For the rest 
of the parameters the Shadow algorithm was set as above. The obtained 
summary statistics show that if the values of A are too high, such as in 
A = (0.1,0.1), then the approximate posterior is rather far away from the 
true posterior. On the other way around, if A = (0.001,0.005) the summary 
statistics seem to indicate a good behaviour of the algorithm. Nevertheless, 
a visual investigation of the sample path confirmed as expected, that in this 
case, the simulated chain is slow mixing. Clearly, there is a compromise in 
choosing the A parameter. A too high value of A will make the algorithm 
diverge, since the distance between the approximate and the true posterior 
linearly depends on A. An algorithm using a too small value of this pa¬ 
rameter will produce samples that indeed may tend to approach slowly the 
true posterior distribution, but the mixing properties of the simulated chain 
are not guaranteed. The alternative intermediate values we have used for 
A show good values for the summary statistics, a proper mixing behaviour 
and independence of the initial conditions Oq. 


4.2 Posterior approximation for the analysis of spatial pat¬ 
terns 

The Shadow algorithm is applied here to the statistical analysis of three 
spatial data sets. The first two data sets consist of several patterns which 
are realisations of a Strauss model and a Candy model, respectively. The 
third one is a a single point pattern which is real cosmological data m- 


4.2.1 Strauss model patterns 

The Strauss model simulates random patterns made of repulsive 

points. The probability density of the process is 

p{y\6) oc = exp(n(y) log /? -t Sriy) log 7 ). (13) 

Here y is a point pattern in the window W, while t{y) = {n{y), Sr{y)) and 
9 = (log/3, log 7 ) are the sufficient statistics and the model parameters vec¬ 
tors, respectively. The sufficient statistics components n(y) and Sr{y) rep¬ 
resent each, the number of points in W and the number of pairs of points 
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MH Sampler: sampling path 



1.5 1.6 1.7 1.8 1.9 2.0 

Theta 1 


MH Sampler - first parameter 



IT) 


0 200 400 600 800 1000 

Sampling time 


MH Sampier - second parameter 



o 

cd 


0 200 400 600 800 1000 

Sampling time 


ABC Aigorithm: sampiing path 


o 

o 


to 

o> 

(N 

nj 

oj o 

E 


to 

cd 


o 

cd 


1.5 1.6 1.7 1.8 1.9 2.0 2.1 

Theta 1 


ABC Aigorithm - first parameter 




1-1-1-1-1-r-* 

0 200 400 600 800 1000 

Sampling time 


ABC Aigorithm - second parameter 



*1 - 1 -:- 1 - 1 - 

0 200 400 600 800 1000 

Sampling time 


Figure 2: Sample path for the Normal posterior. Left colum: the MH 
algorithm results - from top to bottom the joint parameter path, the 9i time 
series and the 02 time series. Right column: the ABC Shadow procedure - 
from top to bottom the joint parameter path, the 0i time series and the 02 
time series. 
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Parameters 

Summary statistics for Normal posterior sampling 
01,02 

0 

<1 

Qb 

Q25 

Qso 

0 

Q75 

Q95 

A = (0.001,0.005) 

1.68 

1.76 

1.82 

1.82 

1.87 

1.95 

00 = (2,9) 

8.41 

8.81 

9.13 

9.09 

9.38 

9.83 

A = (0.I,0.I) 

-0.45 

0.21 

1.96 

1.79 

3.30 

3.94 

00 = (2,9) 

5.87 

6.57 

7.03 

7.06 

7.52 

8.31 

A = (0.01,0.05) 

1.58 

1.69 

1.76 

1.76 

1.83 

1.92 

00 = (2,9) 

8.42 

8.77 

9.03 

9.07 

9.33 

9.82 

A = (0.005,0.025) 

I. 6 I 

1.70 

1.76 

1.76 

1.83 

1.93 

00 = (2,9) 

8.44 

8.78 

9.05 

9.07 

9.32 

9.74 

A = (0.005,0.025) 

1.60 

1.70 

1.76 

1.76 

1.82 

I.9I 

00 = ( 10 , 20 ) 

8.35 

8.78 

9.03 

9.06 

9.33 

9.83 

A = (0.005,0.025) 

1.61 

1.70 

1.77 

1.76 

1.83 

1.92 

00 = (- 10 , 1 ) 

8.39 

8.80 

9.06 

9.08 

9.36 

9.80 


Table 2: ABC Shadow behaviour depending on the A parameter and Oq 
initial conditions. 


at a distance closer than r. The {3 parameter controls the number of points 
in a configuration, while the 7 parameter controls the relative position of 
the points in a configuration. The model is well defined for /3 > 0 and 
7 €]0,1]. If 7 = 1 then the point process is a Poisson point process, since no 
interaction between points exists. If 7 G (0, 1) the points in a configuration 
exhibit a repulsion interaction. The range parameter r is usually considered 
known. It cannot be taken directly into account in our framework since the 
likelihood p{y\0) is not a differentiable function of it. 


Here, it is not possible to make a direct comparison that certifies the quality 
of the results given by the ABC Shadow algorithm. Nevertheless, for some 
particular situations, it may be investigated whether the inference from the 
approximated posterior distribution is close to the inference performed using 
classical methods. 


Let us assume that a random sample of size m made of independent re¬ 
alisations of a Strauss model with parameter 6 is observed. It is easy to 
check following [laila], that the maximum likelihood estimate 9 satisfies 
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the equation 

m 

^ t{yi) - mE^t(X) = 0, 

i=l 

with yi,..., Ym the m point patterns forming the sample. The sum and the 
expectation are computed componentwise. Clearly, the maximum likelihood 
estimate approaches the true model parameters whenever the sample size m 
increases. Since the considered parameter domain is a compact set 0 C 
the likelihood based inference and the one based on a posterior distribution 
build with uniform priors on 0, are equivalent. 

This property can be used to test the ABC Shadow algorithm: within this 
context, the maximum of the approximate posterior should be close to the 
maximum likelihood estimate and the true model parameters. For saving 
computational time, instead of observing m independent samples, an average 
pattern with sufficient statistics given by ^ ^(y*) can be considered. 

In the following, the previous ideas are used to test our method. For this 
purpose, the Strauss model on the unit square W = [0,1]^ and with density 
parameters 13 = 100, 7 = 0.2 and r = 0.1, was considered. This gives for 
the parameter vector of the exponential model 9 = (4.60, —1.60). The exact 
algorithm CFTP (see Chapter 11 in [19]) was used to get 1000 samples from 
the model and to compute the empirical means of the sufficient statistics 
^(y) = (h(y), sF(y)) = (34.33,5.31). The ABC Shadow algorithm was run 
using this sufficient statistics vector as observed data. The prior density 
p(9) was the uniform distribution on the interval [3.5, 5.5] x [—5,0]. Each 
time, the auxiliary variable was sampled using 100 steps of a MH dynam¬ 
ics [l2l[T9]. The A and n parameters were set to (0.01,0.01) and (200,200). 
The algorithm was run for 10® iterations. Samples were kept every 10® steps. 
This gave a total of 1000 samples. 

Figure [3] shows the histograms obtained using the ABC algorithm for sam¬ 
pling the posterior density given by the Strauss model (fTKI) . The kernel 
density estimates are superimposed on these histograms. The MAP (Maxi¬ 
mum a posteriori) is computed by taking the maximum of these estimated 
densities. The obtained value is 9 = (4.63,-1.53). The obtained median 
and mean posterior estimates were (4.606,-1.669) and (4.603,-1.700), re¬ 
spectively. All these values are close to the true values of the parameters. 
Note that in this case, the true maximum likelihood estimates are also equal 
to the true model parameters. Nevertheless, nothing can be said concerning 
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the general shape of the distribution. 



log(beta) log(gamma) 

Figure 3: Strauss posterior sampling using the ABC procedure: histograms 
and kernel density estimates of the marginals. 


4.2.2 Candy model patterns 

The Candy model is a marked point process that simulates connected seg¬ 
ments m- The model was applied with success in image analysis and 
cosmology [ 201122 ] • 

Here, a segment or a marked point y = is given by its centre w, its 

orientation ^ and its length 1. The orientation mark is a uniform random 
variable on M = [0, vr), while the length is a fixed value. The Candy model 
we consider is given by the following probability density 

^(yl^) oc expiOdUdiy) + OgUsiy) + 6*/n/(y) + 9rnr{y)) (14) 

with e = {9d,9s,9f,9r) and t{y) = {ndiy),ns{y),nf{y),nr{y)) the parame¬ 
ters and sufficient statistics vectors, respectively. The parameter 9d controls 
the number nd of segments connected at both of its extremities or doubly 
connected, the parameter 9s controls the number Ug of segments connected 
at only one of its extremities or singly connected and the parameter 9f 
controls the number Uf of segments that are not connected or free. The 
connectivity interaction of two segments is based on the relative position of 
their extremities and also on their relative orientation. Two segments with 
only one pair of extremities situated within connection range Tc and with 
absolute orientation difference lower than a curvature parameter Tc are con¬ 
nected. The parameter 9^ controls the number of pairs of segments that 
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are too close and not orthogonal. This interaction controls the segments ten¬ 
dency to form clusters or to overlap. The orthogonality is controlled through 
a curvature parameter r^. For more details and a complete description of 
the Candy model interactions and properties we recommend m- 

Figure m shows a realisation of the Candy model on W = [0, 3] x [0,1]. The 
segment length is I = 0.12, the connection range is rc = 0.01, and the cur¬ 
vature parameters are Tc = Tr = 0.5 radians. The model parameters are 
dd = 10, 9s = 7, 9f = 3 and 9r = —1. It can be observed that with these 
parameters a rather connected pattern of segments is formed. 



Figure 4: Realisation of the Candy model. 

The same strategy as for the Strauss model was used to test the ABC Shadow 
algorithm. An adapted MH algorithm m was used to obtain 1000 samples 
of the previous model and to compute the vector of the empirical means 
of the sufficient statistics t(y) = (hrf(y) = 51.10,hs(y) = 101.06, h/(y) = 
19.97, hr.(y) = 72.89). These statistics were considered as the observed data 
for our ABC Shadow algorithm. The prior density p{9) was the uniform 
distribution on the interval [2,12] x [2,12] x [2,12] x [—7,0]. Each time, 
the auxiliary variable was sampled using 2000 steps of an adapted MH dy¬ 
namics. The A and n parameters were set to (0.01,0.01,0.01,0.01) and 
(500, 500, 500, 500). The algorithm was run for 10® iterations. Samples were 
kept every 10^ steps. This gave a total of 1000 samples. 

Figure [5] shows the histograms obtained using the ABC algorithm for sam¬ 
pling the posterior density given by the Candy model (I14p . The kernel 
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densities estimates are superimposed on these histograms. The MAP is 
computed by taking the maximum of these estimated densities. The ob¬ 
tained value is 9 = (9.96, 6.99, 3.02, —1.00). The obtained median and mean 
estimates were (9.995, 7.005, 2.977, -1.014) and (9.998, 7.008, 2.975, -1.018), 
respectively. Again, the obtained values are close to the true parameter val¬ 
ues. 



theta_d 


thetaj 



theta_s 


lheta_r 


Figure 5: Candy posterior sampling using the ABC procedure: histograms 
and kernel density estimates of the marginals. 


4.2.3 Real data application 

The galaxies in our Universe are not uniformly distributed. They tend to 
form structures such as filaments, sheets and clusters [T6] . Figured shows 
a sub-region from a two-dimensional cosmological catalogue [22]. The data 
set contains a number of 163 points that represent the galaxy positions in 
the considered region. The observation domain is scaled to a square window 
W = [0,1] X [0,1]. The visual inspection of the data set indicates a cluster- 
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ing tendency of the pattern. In the following, it will be shown how the ABC 
Shadow algorithm can be used to test this hypothesis. 



Figure 6: Sample of a cosmological data set. The points represent galaxy 
positions in a region of our Universe. 

Such point patterns are typically analysed using summary statistics [HKH]. 
For a stationary point process Y of intensity A > 0, some of the most used 
summary statistics are: 

• the K function, which is proportional to the expectation under the 
reduced Palm distribution of N{b{o, u)), the number of points in a ball 
of radius u centred at the origin o, that is a point of the considered 
stationary process which is not counted 


XK{u) =E|, [N{b{o,u))] 
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• the nearest neighbour distance distribution function 

G{u) =F[{N{b{o,u) > 0)) 
with Pq the reduced Palm distribution, 

• the empty space function 

F{u) =P(Yn6 (o,u) / 0) 


with P the distribution of Y, 


• the J function which compares nearest neighbour to empty distances 


J{u) 


1 - G{u) 
1 - F{u) 


for all M > 0 such that F{u) < 1. 


These summary statistics have exact formulas for the stationary Poisson 
point process with intensity A ; 


K{u) = TTU^, 

F{u) = 1 — exp[— 
G{u) = F(u), 

J{u) = 1. 


So, for a given point pattern, the estimators of these statistics give indica¬ 
tions about how far from the Poisson point process, the considered pattern 
is. For the K and G statistics, higher values of the observed summaries than 
the theoretical Poisson statistics suggest clustering of the pattern, while the 
lower ones suggest repulsion. For the F and J statistics, it is the other way 
around: lower observed values than the theoretical Poisson ones indicate 
clustering, while the higher values recommend repulsion. To better answer 
this type of questions for an observed point pattern, an envelope test can be 
done. This test simulates a Poisson point process with intensity estimated 
from the observed pattern, and for each simulated pattern the summary 
statistics are computed. At the end of the simulation, Monte Carlo confi¬ 
dence envelopes are created. So, if the observed summary statistics belong 
to the envelope, the null hypothesis that the observed point pattern is the 
realization of a Poisson process is not rejected, hence the point pattern ex¬ 
hibits a completely random structure. For a thorough presentation of the 
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point process summary statistics and of the envelope tests based on them, 
we recommend muH] and the references therein. 

Figure [7| shows the result of the envelope for the galaxy point pattern. All 
the statistics exhibit the clustered character of the pattern. Despite the fact 
this very useful analysis result may be considered reliable, there is no in¬ 
formation concerning the “strength” of the clustering character. This may 
be achieved by assuming the pattern as the realisation of known model and 
estimate its parameters. Still, parametric inference using these summary 
statistics should be done with care, since it may be misleading. For in¬ 
stance, the authors in proved the existence of point processes that are 
not Poisson, but exhibiting K{u) = ttu^ or J{u) = 1. Another possibil¬ 
ity to be considered is to perform pseudo-likelihood (PL) or Monte Carlo 
maximum likelihood (MCML) estimation. Now, the PL estimation lacks of 
precision if the objects interaction is too strong, while the MCML compu¬ 
tation requires re-sampling of the model whenever the initial conditions are 
too far away from the true maximum likelihood. As an alternative to these, 
we use posterior sampling using the ABC Shadow algorithm. 

The area-interaction process was introduced by [1] in order to produce clus¬ 
ter point patterns. The probability density of the process is 

p{y\e) oc = exp(re(y)log/ 3 -kar(y)log 7 ). (15) 

Here y = {wi,W 2 ,..., is a point pattern in the finite window IT, while 
t{y) = (’^(y)) «r(y)) and 6 = (log/I, log 7 ) are the sufficient statistics and 
the model parameters vectors, respectively. The statistic re(y) represents 
the number of points in a configuration. The statistic ar{y) is given by 

^ vrr^ vrr^ 

and it is proportional to the surface of the union disks of radius r centred in 
the points of y. It can be interpreted also as the number of disks needed to 
“cover” the area of the disk pattern Ar{y) = ^2=ib{wi,r). Here, as outlined 
by [4|, there is a direct link between the statistic ar{y) and the ’’reduced 
sample” estimate of the empty space function F(r), in the sense that the 
vector (n(y),F(r)) is also a sufficient statistics vector. As for the Strauss 
model, the (3 parameter controls the number of points in a configuration, 
while the 7 parameter controls the relative position of the points in a con¬ 
figuration. The model is well defined for /3 ,7 > 0. If 7 = 1 then the point 
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F envelope test G envelope test 



Figure 7: Summary statistics envelope test (100 simulations) for the galaxy 
pattern: shaded (gray) region - the Monte Carlo envelopes, dotted line (red) 
- the theoretical statistics, continuous (black) line - the observed statistics. 
These results (estimators computations and simulations) were obtained us¬ 
ing the R spatstat library [3]. 
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process is a Poisson point process, and no interaction between points exists. 
If 7 > 1 then the configurations inducing disk patterns with large areas will 
be penalised, hence producing a clustering effect of the points in the pattern. 
If 7 < 1, then the configurations inducing disk patterns with small areas will 
be penalised, hence generating a repulsion effect. As for the Strauss process, 
the range parameter r is considered known, and it cannot be taken directly 
into account because the likelihood p{y\0) is not a differentiable function of 
it. 

In the following, we explore the structure of the galaxy pattern by assum¬ 
ing it is a realisation of an area-interaction process and by investigating its 
posterior distribution. Within this context, the posterior distribution gives 
important information related to the direction and the strength of points 
aggregation in the pattern. If the prior used is the uniform distribution, as 
in the previous situations, the maximum of the posterior is the maximum 
likelihood estimate. In this case, the posterior distribution gives, simultane¬ 
ously, information about the general morphology of the pattern and about 
the model able to reproduce a pattern having sufficient statistics similar to 
the observed ones. 

For this purpose, the ABC Shadow algorithm was used in the following way. 
First, similarly to the summary statistics analysis, a finite set of values for 
the range parameters was fixed. From the galaxy point pattern, for each 
fixed value of r, the sufficient statistics of the area-interaction model were 
computed. Table [3] presents these values. 


Data for the Galaxy pattern 

r 

0.01 

0.02 

0.03 

0.04 

0.05 

0.06 

0.07 




n{y) = 

163 




-ar{y) 

135.91 

114.05 

96.44 

82.23 

69.85 

59.05 

49.73 


Table 3: The observed sufficient statistics computed for the galaxy pattern, 
depending on the range parameter r. For all these parameters n(y) remains 
constant, while ar{y) depends on r. 

For each r value an ABC Shadow algorithm was initialised with the com¬ 
puted sufficient statistics as the corresponding observed data. Except, the 
observed sufficient statistics, the parameters of all the algorithms were all 
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the same. The prior density p{9) was the uniform distribution on the inter¬ 
val [2,12] X [—5, 5]. Each time, the auxiliary variable was sampled using 250 
steps of a MH dynamics. The A and n parameters were set to (0.01,0.01) 
and (100,100). The algorithm was run for 10® iterations. Samples were kept 
every 10® steps. This gave a total of 1000 samples. 

Figure [ 8 ] shows the boxplots of the posterior distributions for the parameters 
log (5 and log 7 conditionally on their corresponding range parameter. Since, 
the boxplots whiskers for the logy posterior are far away from 0 it may be 
concluded that there is a strong evidence of clustering of the pattern, for 
the considered ranges. 


Posterior of log (beta) parameter 


Posterior of log (gamma) parameter 



Range r = 0.01,0.07 


Range r = 0.01,..., 0.07 


Figure 8 : Box-plots of the posterior distributions for the parameters of the 
area-interaction process estimated from the galaxy pattern, given different 
values for the interaction radius. 

The model parameters can be estimated for each radius. As in the previous 
examples, the MAP estimate equals the Maximum Likelihood (ML) esti¬ 
mate. In order to verify these results, the estimation errors are computed as 
in m- The results are shown in Table 01 For their computation, a simula¬ 
tion of 10000 samples for each estimated model was performed. The asymp¬ 
totic standard deviation approximates the difference between the maximum 
likelihood and the true model parameters, which are both unknown. The 
Monte Carlo standard deviation indicates here the difference between the 
MAP estimate and the maximum likelihood. 
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Estimation errors 

r 

0.01 

0.02 

0.03 

0.04 

0.05 

0.06 

0.07 

Asymptotic standard deviation 

O'log^ 

0.20 

0.17 

0.13 

0.11 

0.10 

0.10 

0.08 

^log7 

0.26 

0.28 

0.27 

0.30 

0.34 

0.40 

0.52 

1 

Vlonte Carlo standard c 

leviation 

^log § 

0.001 

0.002 

0.001 

0.002 

0.002 

0.002 

0.003 


0.002 

0.003 

0.004 

0.006 

0.008 

0.012 

0.024 


Table 4: Estimation errors computed for the MAP estimates obtained from 
the galaxy pattern. 


5 Conclusion and perspectives 


The proposed method enlarges the panel of existing ABC techniques. To our 
best knowledge, it is the first time an ABC based statistical analysis is based 
on these point processes. The method was tested and tuned on known pos¬ 
teriors from the exponential family models. The results obtained on posteri¬ 
ors of point processes are coherent with the summary statistics analysis and 
the maximum likelihood estimation. Compared with the pseudo-likelihood 
based inference, our algorithm has the advantage that its output is a dis¬ 
tribution that tends to approach the true posterior. The theoretical results 
hold for a class of models larger than the exponential family and they allow 
the use of different prior distributions. Clearly, the choice of the appropri¬ 
ate statistics and the construction of the test procedures for the algorithm 
within this general context, are open problems. 


As for all the other methods of sampling posteriors, the simulation of the 
auxiliary variable should be done exactly. Nevertheless, the numerical re¬ 
sults obtained with less perfect samplers were satisfactory. The strong point 
of this method is that by getting close to the posterior distribution, no ’’use¬ 
less” samples are produced. As perspectives, we mention range parameters 
estimation and model validation. 
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Appendix : proof of the results 

Proof of Theorem [1] 

(i) Both density functions vanish outside the ball b{9, A/2). For ijj G 
6(0, A/2), the integral mean value theorem applied to the denominator 
of (7A leads to 


/(x|b) 



for some 9* G 6(0, A/2). The positivity and the continuity of the 
density p (uniform continuity indeed, since 0 is compact) allows us to 
do the following. Let m(x) := inf,^g©p(x|0) > 0 since it is actually a 
minimum. For A G T© we have 



A 


where the last supremum is independent of ijj and 9*, and approaches 
to 0 as far as A does. With the regularity condition on p(x|-), we can 
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tune up the inequality, using the (differential) mean value theorem, to 


\qa{G 'ij^) — Ua{ 6 dip < m(x) sup \\DQp{x.\'ij;* 

i/i*G0 


:= C'i(x,p,0)A 

where Ci(x,p, 0) is a constant depending on x,p and 0. 

(ii) As previously, the use of the integral mean value theorem gives the 
result, since 


sup 

■0G0 


qA{9 V^|x) 
QAiip ^Ix) 


^V(»,A/2)(V) 


sup 

ip&{e,A/‘i) 


/(x|e*) N 

c(e*) ) 

/(x|V>*) 2, 


/(x|V’) 

c(^) 

/(x|e) 

<e) 


< 


sup 

V>Gfc{6»,A/2) 


/(x|V’) 

/(x|e) 

'wr 


fwr) 

c[ip*) 

/(x|g*) 

-T(0^ 


1 < 


M(x)m(x) ^ sup 


/(x|V^*) 

c(V'*) 


/(x|r) 

0(0*) 


where M(x) := sup 0 g 0 p(x| 0 ) < oo is a maximum and 0* € 6(0, A/2), 
pj* G A/2) are obtained from the integral mean value theorem. 
As in (i), under the regularity condition on p(x|-), this inequality can 
evolve to 


sup 

ip€B 


gA(0 V^|x) 

QAiip 0|x) 




< M(x)m(x) ^A sup ||D 0 p(x|V^*)|| 
1/1*60 


:= C' 2 (x,p, 0)A 


where C 2 (x,p, 0 ) is a constant depending on x,p and 0 . 


Proof of Proposition [T] 

If n = 1, the dehnition of the transition kernels, the introduction of the 
term Ua{ 9 'ip)oii{6 ip) — Ua{6 —>■ 'ip)ai{9 —>■ '*/’)) then the use of the 
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triangle’s inequality and also the boundedness of functions “*(') 

as(-), allow us to write 


\P,{d,A)-Ps{d,A)\ < 

\qA{0 i’)oii{9 Ip) — Ua{ 9 'p)as{0 -)■ 'p)\d'ip + 


'V’66(6l,A/2) 


+ 

< 3 
+ 2 


u{e) [ 

J 'll) 


\qA{d - ai{9 V’)] - UAi9 -)• - as{9 


V'6fe{6»,A/2) 

\qA{9 ip) - Ua{ 9 'ip)\d'ip+ 


/i/'eb(0,A/2) 


/ 

J 'll) 


i/>eb(0,A/2) 


[/ a ( 6 ' —>■ 'ip)\as{9 ip) - as{9 -)> ip)\dip 


Under the hypothesis of Theorem [U and then applying Theorem [T]^i) and 
Corollary [H the transition kernels of the ideal and the shadow Markov 
chains, respectively are uniformly close as well (say for any e > 0 , there 
exists A(e, 1) > 0 such that we have \Ps{9,A) — Pi{9,A)\ < e if A < A(e, 1), 
independently of 9 and A). 

If p(x|-) G C^(0), then the previous inequalities can be completed to 


\Pi{9,A) - Ps{9,A)\ < 3C'i(x,p, 0)A + 2 C' 2 (x,p, 0)A := C' 3 (x,p, 0)A 

giving a candidate expression for Ao(e, 1) := 6 * 4 (x,p, 0)e. The constants C 3 
and 6*4 depend on x,p and 0 . 


For n > 1 we get by induction that 

\p(n) _ pjn)| < |pW _ p,(^-i)pji)| p ip.l’^-bpji) _ p0)| 

< - Pi’""^)| 

< ... < _pa)| p g 

uniformly in 9 and A for all A such that 


A < Ao(e,n) := Ao(e/n, 1) = C' 4 (x,p, 0)-. 

n 


( 16 ) 
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